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CHAPTER 1 


INTRODUCTION 

1.0 Background 

The NASA Langley Research Center is engaged in obtaining in- 
flight nearby and direct strike lightning data on an instrumented F-106B 
thunderstorm research aircraft. The aircraft is equipped with sensors 
which measure surface electric fields, surface magnetic fields, and 
current on the pitot boom with a data system incorporating a transient 
digitizer. The objective of this research is to characterize the 
lightning environment in such a way that it is applicable to other air- 
craft as well. The aircraft has encountered numerous direct lightning 
strikes during the summers of 1980, 1981, and 1982, and has likewise 
recorded numerous electromagnetic transient waveforms. 

In order to derive the lightning environment from these measure- 
ments, the effects of the aircraft have to be removed from the data to 
yield the environment which would be there without the aircraft present. 

The objective of the research discussed in this report is to develop a 
methodology to determine this environment and to apply it to the existing 
data. This research is a continuation of a previous effort reported in 
[1]. in which an initial data interpretation approach was developed. 

This reference provides much of the background for the results pre- 
sented here. This includes a discussion of the three dimensional 
finite difference technique used to model the electromagnetic response 
of the F-106B, preliminary work in nonlinear analysis, and an initial 
interpretation methodology development. 

The following general research areas are discussed in this report: 

1. Review of 1981 and 1982 Data (Chapter 2) 

2. Interpretation Methodology Development (Chapter 3). 

3. Nonlinear Air Breakdown Modelling (Chapter 4). 

4. Model Validation Studies (Appendix A). 
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Item 4 above was published as a conference paper and is included 
in Appendix A as such. Finally, a discussion of the overall results and 
suggestions for improving the test program are given in Chapter 5. 

1 . 1 Summary 

The emphasis in this research has been the development of method- 
ologies and tools which can be applied to the data interpretation problem. 

One of these tools is a linear approach for obtaining the lightning channel 
current in the absence of the aircraft as a function of various channel para- 
meters. When applied to only one measurement point, however, the result is 
not unique. Initial attempts to obtain a unique answer from correlated 
1982 measurements met with limited success, indicating that the proper 
channel parameters were not included. 

The concept of finding ratios of Fourier transforms of simul- 
taneous responses was introduced with the objective of overcoming the 
problem concerning uniqueness. Although this concept needs to be further 
developed and applied, initial results are encouraging. 

It should be emphasized that it is not certain that linear 
analysis is sufficient to properly understand the in-flight data. With 
this in mind, a nonlinear air breakdown model based on first principles, 
was developed and applied. The results indicate that it is entirely 
possible that most of the data obtained thus far can be interpreted as 
the nonlinear attachment of a leader channel to the aircraft. 

Future efforts will concentrate on application of these tools 
to the measured data, and further exploitation of simultaneous measure- 
ments taken in 1982. 
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CHAPTER 2 


REVIEW OF INFLIGHT DIRECT STRIKE DATA 

2.0 Background 

In this chapter, a review of direct strike inflight data from 

1981 and selected data from 1982 is given. This includes Fourier trans- 
form analysis to identify resonances and nulls. Simultaneous data from 

1982 is discussed along with its usefulness in helping to solve the problem 
of uniqueness in data interpretation. 

2.1 Fourier Transforms of Recorded Data 

Fourier transforms are presented for data gathered both during 

1981 and 1982. Transforms are presented in a log-linear format in Figures 

2.1 - 2.17. The vertical axis is dB, where dB = 20 log^g | F((d) |, where 
F(w) is the Fourier transform expressed in MKS units of B(t) or D(t), 

as appropriate. The Fourier transforms were calculated using the algorithm 
described in Appendix B. The horizontal scale is in frequency from zero 
to 30 MHz. Although the measured data has a useful bandwidth up to 

50 MHz, which was used for the analysis done in Chapter 3, it was decided 
to present the data here with a smaller bandwidth and on a linear scale 
because it was found that this provided the best resolution of the aircraft 
resonances . 

In the interest of brevity, the time domain data is not presented 
here. The 1981 data can be found in a published report [2], and the 

1982 data will be contained in a future report. 

The captions require some explaining. The first 7 characters 
give the flight and run number, and the 8th digit indicates whether the 
measurement is B or D. For example, 81-042R1B means that the data shown 
is B for flight number 42, run 1, during the 1981 thunderstorm season. 

Data for 1981 is given in Figures 2.1 - 2.5, and selected data 
for 1982 is given in Figures 2.6 - 2.17. The latter data is important 
because it consists of correlated B and D records for the same lightning 
event, a feature which will be discussed later in more detail. 
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Figure 2.3 81-043 R2B Figure 2.4 81-043 R3B 
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Figure 2.17 82-039 R2D 


Tables 2.1 and 2.2 summarize the dynamic range, the inverse 
(time domain) record length, and the exact frequencies of the nulls and 
peaks for the 1981 and 1982 data, respectively. The first column describes 
the event, and the second column gives the estimated dynamic range, which 
is 20 log^g(N), where N is the maximum number of discrete levels used in 
the time domain data record. This number is an indication of the dynamic 
range of the Fourier transform data as well. It is observed that the 1982 
data has a somewhat better dynamic range. 

The nulls in the data, especially the ones at the lower 
frequencies, are of interest because if they can be shown to occur 
because of aircraft resonance effects, something might be learned 
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TABLE 2.1 

1981 B Data Summary for Direct Strikes 



EVENT 


DYNAMIC 

RANGE 

dB 

[RECORD . 
LENGTH]' 1 
MHz 


NULLS 

MHz 



PEAKS 

MHz 


81 

-042 

R1 

B 

14 

4.3 

2.85 

. 14.9 


6.25, 

, 11.4 

, 20.9 

81 

-043 

R1 

B 

14 

4 

2.77 

, 9.75, 

14.6 

6.65, 

11.4 

, 17.3 

81- 

-043 

R2 

B 

17 

5.9 

.86, 

15.4 


12.5, 

20.4 


81- 

-043 

R3 

B 

19 

2 

2.0, 

9.3, 19 

.2 

6.4, 

11.2, 

21.7 

81- 

-043 

R4 

B 

17 

5.3 

3.0, 

15 


11.9, 

19.6 



about the channel impedances and attach point locations. For example, 
identification of X/4 resonances could give indications that the air- 
craft has either high or low impedance channel loading, depending 

• • 

upon whether the resonance is observed in D or B, and the attach point 
locations. It is of great interest to determine if the nulls are 
caused by aircraft resonances or if they are related to the record 
length. The inverse of the time domain record length was computed for 
correlation with these nulls. It is observed that there is no obvious 
correlation of the lowest frequency nulls with the record length, 
although this does not prove that they are not related. For example, 
consider the two Fourier transform pairs in Figure 2.18. Clearly, the 
nulls in the amplitude spectrum occur at frequencies related differ- 
ently to the pulse width. Other transforms were computed for various 
arbitrary triangles and other simple waveforms, which showed that the 
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TABLE 2.2 

1982 B and D Data Summary for Direct Strikes 


EVENT DYNAMIC [RECORD NULLS PEAKS 

RANGE LENGTH]' 1 MHz MHz 

dB MHz 


82-037 R4 

D 

25 

1.37 

1.33, 2.86, 6.4, 
8.99 

1.88, 3.39, 5.5, 
7.5, 10.1, 17.8, 
20.0 


B 

24 

1.5 

1.43, 3.75, 6.33, 
16.4 

2.54, 5.19, 7.6, 
10.2, 12.4, 17.8, 
20.2 

82-038 R2 

D 

34 

.71 

2.85, 11.4, 14.9 

6.0, 12.0, 16.2 


B 

30 

1.06 

1.30, 10.1, 14.9 

5.9, 16.0, 20.2 

82-038 R4 

D 

23 

1.1 

1.22, 2.77, 3.87, 
9.4 

2.17, 3.3, 7.54 


B 

23 

1.9 

1.70, 5.46, 9.23, 
16.4 

7.27, 11.1, 18.3 

82-038-R7 

D 

29 

1.2 

1.27, 4.25, 10.8 

1.88, 6.1, 20.2 


B 

22 

1.37 

3.1, 5.0, 7.0, 
9.0, 11.1, 13.1, 
15.0 

2.0, 4.1, 6.1, 

8.0, 10.0, 12.1 

82-038-R8 

D 

21 

5.3 

3.84, 10.45 

7.2, 13.6, 20.0 


B 

21 

4.8 

15.6 

7.4, 19.4 

82-039 R2 

D 

20 

5.0 

3.7, 10.4, 16.8 

7.1, 13.9, 19.9 


B 

26 

4.2 

2.83, 9.75, 15.5 

6.83, 12.1, 20.5 
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Figure 2.18 Rectangular Pulse and Isosceles 
Triangle Transform Pairs 


resulting nulls cannot be simply related to the record length, but 
are related to other peculiarities of the waveforms. Thus it is be- 
lieved that the lower frequency nulls in the data cannot be proven 
to be related to aircraft resonances, and must be assumed to be 
caused by some other source. In addition, the lower frequency nulls 
do not correlate with the basic X/4 resonance at about 3.5 MHz. 

The peaks in the data are shown in the last column. Because 
of the complex nature of the Fourier transform waveforms, it is 
difficult to tell whether a peak occurs because of an aircraft re- 
sonance, or simply because it is between two nulls. The converse is 
also true. In spite of this however, there are certain frequencies 
which consistently appear in much of the data. 
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The 1981 B data show fairly consistent peaks at 6-7 MHz, 

11-12 MHz, and 20-21 MHz. This is in general agreement with the 
resonances reported by Trost [3] for the 1980 data. There is also 
a rather consistent null at 'vlO and ^15 MHz. 

The 1982 data is similar in character. Peaks in the vicinity 
of %7, 12 and ^20 MHz are common. Nulls at ^10 and %15 MHz are also 
evident. 

It is likely that these resonances are related to aircraft 
dimensions. For example, the half wavelength resonance based on 
structural length of the fuselage including the pitot boom is 7.0 MHz, 
and is 8.8 MHz without the pitot boom. The wing tip to wing tip 
resonance is on the order of 13 MHz. One would expect that the fatness 
of the aircraft would reduce these frequencies on the order of 5-10%. 

The null data and peak data shown in Tables 2.1 and 2.2 is 
based on a Fourier transform technique which assumes a straight line 
interpolation between measured data points. The frequencies indicated 
are accurately determined from the computer numerical output directly. 

The frequency of the deepest null or largest peak is underlined. In 
the case in which there is more than one frequency underlined, then 
the depths or peaks at these frequencies are within 10%. The nulls and 
peaks are given only for frequencies less than 20 MHz, in order to 
keep the data set small and to emphasize the primary and first few 
resonances. It should be noted that not all peaks and nulls are given, 
but only those that appear to be significant, which involves some 
judgment. 

Inspection of the transforms in Figures 2.1 - 2.17 results 
in some interesting observations. First, there is a great deal of 
variety in the transform structure. The data for 1981 (Figures 2.1 -2.5), 
for example, has a few broad peaks and valleys, whereas the 1982 data 
(Figures 2.6- 2.17) shows much more structure with considerably more 
peaks and nulls. This difference may be caused by the lower dynamic 
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range of the 1981 data with respect to the 1982 data. Another 
possible cause might be that the 1982 strikes were encountered 
at much higher altitudes ( £ 9 km) than were the 1981 strikes 
(%6 km). For example, 82-037R4D and 82-038R4D show periodic 
nulls and dips every 2.3 - 2.6 MHz. Similarly, 82-038R4D and 
82-038R4B show this pattern at 3.2 - 3.6 MHz intervals. Again, 
the same type of pattern is observed in 82-038R7D and 32-038R7B 
at 2 MHz intervals. This latter feature is probably caused by 
the presence of a second pulse in the time domain waveforms 
delayed 490 ns from the initial pulse. These responses can be 
contrasted with the much less complex structure of 82-038R8D and 
82-038R8B. 

2.2 Correlated 1982 Data and the Problem of Uniqueness 

The data shown for 1982 is of special interest because it 
involves simultaneous B and D records. These records are not corre- 
lated accurately enough to identify which one began first, but 
they are correlated to the extent that they can be identified as 
being caused by the same lightning strike event. Because whatever 
lightning environment is postulated must simultaneously cause these 
two responses, it is expected that confidence in the uniqueness of 
that environment is much greater than that for one derived from a 
single measurement point. 

For a particular channel configuration (impedance, velocity 
of propagation, attach points, etc.) transfer functions Tq (w) and 
Tg (w) can be defined as 


and 



( 2 . 1 ) 


Tg(co) 



( 2 . 2 ) 
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where I (u>) is the channel current, and linearity is assumed. It is 
easy to show then that 


DU) = T D (w) 
B(o>) T B ( “^ 


(2.3) 


If D(w) and B(a )) are measured quantities, and Tp(o>) and TgU) are 
determined from analysis of aircraft response to certain channel configu- 
rations, one should be able to reduce the uncertainty in the uniqueness 
problem by finding the channel configuration given TAU)/T; U) which most 
nearly matches the measured ratio D(w)/B(u>). 


These ratios must be identical at all frequencies, which at 
this point seems to present a difficult problem in view of the possible 
large number of combinations of channel parameters and attach points. 
Nevertheless, this question must be addressed. 

The response ratio RR can be defined in logarithmic 

form as 


RR = 20 log|l^-|> 

D(u>) 


(2.4) 


with D and B expressed in MKS units. 

Figures 2.19 - 2.24 are plots of these ratios determined from 
the measured data. It is noted, of course, that they are not identical, 
although there are certain peak frequencies which are evident in much 
of the data. The differences between the figures arise from different 
channel parameters, attach points, or nonlinear effects. A summary of 
the peaks for the measured waveforms, as well as for some computed 
waveforms, is given in Table 2.3. 

Because the response at low frequencies is easier to under- 
stand than that at high frequencies, it seems reasonable to assume 
that the low frequency response ratio could be indicative of attach 
point locations. That is, one would expect this ratio to be different 
for nose-tail, nose-wing, or wing-wing configurations, for example. 
However, it is found that the dynamic range (or alternately, the 
granularity of the amplitude resolution of the digitized data) is not 
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Figure 2.23 RR for 82-038 R8 Figure 2.24 RR for 82-039 R2 


sufficient for accurate determination of the low frequency levels. 

This is particularly true for B data, which is almost always bipolar 

and has the property that the / Bdt (which is the same as the low 

o 

frequency value of the Fourier transform of B) is very close to zero, 

that is, it has a small DC value. This DC value is of the same size 

(or smaller) as the value one would obtain by assuming an error of 

1/2 of the amplitude digitization resolution over the entire record 

length. The latter value is a measure of the absolute best accuracy 

one could obtain with the discrete data system. Therefore, errors in 

B at low frequencies were computed to be in the range of -°° dB (which 

corresponds to/ Bdt < max digitization error) to + 10 dB. Errors 
o . 

in DC values of D are smaller, because D tends to have a significant 
DC value. The errors for D were computed to be in the range of -4 dB 
to +3 dB. The results are summarized in Table 2.4. The second column 
shows the response ratio at DC as obtained directly from the Fourier 
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TABLE 2.3 


Summary of Peaks for Response Ratios 


EVENT 

FREQUENCY, MHz 

RATIO dB 

82-037 R4 

2.854 

53.2 


4.348 

46.8 


6.382 

60.1 


9.00 

52.5 


11.44 

48.5 

82-038 R2 

1.792 

40.3 


2.845 

58.6 


3.574 

52.4 


5.167 

44.6 


10.00 

43.8 


11.44 

53.7 

82-038 R4 

1.207 

44.9 


2.782 

44.4 


3.871 

64.5 


4.915 

49.0 


9.568 

44.2 


11.17 

44.7 

82-038 R7 

1.27 

52.9 


2.512 

45.0 


4.213 

52.0 


8.713 

44.7 


10.72 

53 

82-038 R8 

3.889 

49.4 


10.45 

44.1 

82-039 R2 

3.727 

63.8 


10.54 

42.5 

Computed linearly: Nose Entry 

2.548 

42.2 

Tail Tip Exit, Z Q = 100 0, 
v p = 3 x 10 8 , R/sl = 0 

12.52 

41.7 

Computed linearly: Nose Entry 

2.305 

50.2 

No Exit, 1 Q = 100 0, 
v p = 3 x 10 8 , R h = 50 n/m 

11.2 

41 .4 
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TABLE 2.3 (CONT'D) 

Summary of Peaks for Response Ratios 


EVENT FREQUENCY, Mhz 

Computed Nonlinearly 

A. Leading Edge of Left Wing 1.594 

I = -100A 3.817 

5.338 
6.103 
7.372 
9.57 

10.5 

11.7 

12.61 

25.03 

Computed Nonlinearly 

B. Leading Edge of Left Wing 1.522 

I = +100A 2.674 

4.339 
6.418 
7.03 
8.578 

10.45 

Computed Nonlinearly 

C. Nose, I = -100A 2.602 

4.411 

6.571 

8.605 

10.99 

12.34 

Computed Nonlinearly 

D. Leading Edge of Right Wing 

I = + 100 A 2.197 

4.312 

6.382 

8.308 

9.89 

11.89 

27.8 

Computed Nonlinearly 

E. Right Side of Tail 1.81 

I = -100A 4.339 

6.877 

9.244 

11.53 

13.06 

26.83 


RATIO, dB 


24.0 

29.9 

35.0 

35.0 
40.4 

40.1 

36.9 

39.2 

42.8 

51.9 


18.6 

30.8 

44.8 
44.8 
38.4 
40.2 
38.33 


44.1 

37.0 
36.3 
36.9 

40.0 
33.8 


40.4 

52.1 

50.2 

45.4 
46.0 

41.5 
65.4 


36.9 

39.9 
41.24 

42.5 

49.5 
59.3 
70.8 
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TABLE 2.4 


Low Frequency Response Ratios for Selected 1982 Simultaneous Data 


EVENT 

RESPONSE RAT 1C 
at DC (dB) 

POSSIBLE VARIATION 
(dB) 

82-037 R4 

28.3 

(-16.9, + 8.5) 

82-038 R2 

19.4 

(-8.9, + 6.2) 

82-038 R4 

21.4 

(-, + 9) 

82-038 R7 

25.0 

(-29, + 8.3) 

82-038 R8 

22.2 

(-», + 14.4) 

82-039 R2 

33.2 

(-8.2, + 5.9) 


transforms. The third column shows the possible variation in this 
ratio, as determined from the error analysis of B and D. The symbol ® 
means that the possible error is larger than the calculated DC value, 
thus allowing the possibility that the DC value is zero. Even apart 
from those two special cases, the allowable spread can be as great 
as 37 dB, which clearly indicates that the data is not sufficiently 
accurate to allow the use of DC levels to infer attach points. This 
problem can be helped by increasing the dynamic range of the in- 
strumentation. It is understood that there are plans to replace the 
current 6 bit digitizer with an 8 bit unit, which will allow up to a 
maximum of 12 dB increase in the dynamic range. 

It thus seems that if the response ratio's DC values cannot 
be used to help in the uniqueness problem, then the peaks of the 
response ratios may provide some insight. The response ratios have been 
calculated for two linear channel configurations indicated in Figures 
2.25 and 2.26 and for five nonlinear attachment cases shown in Figures 
2.27 to 2.31. These nonlinear cases correspond to those of section 
4.6 with 6 % water content and a relative air density of .5. These 
can be compared directly with the measured ratios in Figures 2.19 - 
2.24. Table 2.3 also gives the peak values for the computed ratios. 
Several items are worthy of mention. 
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Linearly Computed RR for Nose 
Entry and No Exit, Z 0 = 100 p, 

R/l = 50 p/m, v = 3 x 10 8 m/sec 
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Frequency MHz 

Figure 2.28 

RR for Nonlinear Attachment to Leading 
Edge of Left Wing I = +100 A 



Frequency Hlz Frequency IHz 


Frequency mz 

Figure 2.27 

RR for Nonlinear Attachment to Leading 
Edge of Left Wing, I = -100 A 



Figure 2.29 Figure 2.30 

RR for Nonlinear Attachment to RR for Attachment to Leading Edge 

Nose, I = -100 A ^ of Right Wing I = +100 A 
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Figure 2.31 

RR for Nonlinear Attachment to 
Right Side of Tail, I = -100 A 


First, the computed ratios have a larger dynamic range than 
the measured ones. This is so, however, because the computed data has 
a larger dynamic range (> 100 dB) than does the measured data (20-30 dB), 
and thus no significance can be attached to this difference. 

Second, the measured ratios often have considerably more 
structure in them than do the linearly computed ratios, and the non- 
linearly computed ratios have more structure in them than do the linearly 
calculated ones. 

Third, the peaks for measured ratios are not free from noise 
problems either. In many cases, the peaks occur because of nulls in 
|D(o))| . The depth of these nulls is greater than what the dynamic range 
of the original data would allow. Thus, the validity of the amplitude 
of the peaks has to be determined on an individual basis by studying 
the original waveforms. 
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Fourth, there are some modest agreements between the 
computed and some of the measured ratios. For example, the ratio for 
82-038R8 has features quite similar to those of Figure 2.26. Also, the 
computed peak frequencies of about 11 and 2.4 MHz are evident in 
several of the measured waveforms. Thus there appears to be some hope 
for using this type of technique to help resolve the uniqueness problem. 

Fifth, there is a great variety in structure of both measured 
and computed ratios. This is particularly true for the ratios obtained 
from responses computed nonli nearly. 

Finally, the utility of the response ratios has not been fully 
evaluated at this point. Work in this area began near the end of this 
research period, and more work is planned for future efforts. In par- 
ticular a whole catalogue of ratios for various attach points and 
channel parameters needs to be developed, and ratios from other 1982 
correlated data sets need to be computed. 
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CHAPTER 3 


INTERPRETATION METHODOLOGY DEVELOPMENT 
3.0 Background 

The response of an aircraft to a lightning strike can be 
significantly affected by the character of the lightning channel. 

Consider as an example two different channels containing a given current, 
one with a high characteristic impedance and the other a low impedance. 
When the high impedance channel interacts with an aircraft, the current 
in the channel is not greatly changed, and the current flowing onto 
the aircraft is approximately the original value in the channel. When 
the low impedance channel interacts, the current in the channel can be 
changed significantly, and the current flowing onto the aircraft is 
usually larger than the initial value in the channel. Since this 
current is the primary source for the aircraft response, it is clear 
that the channel impedance must be considered when studying the 
aircraft response. 

Another characteristic of the lightning channel is the 
propagation velocity of the current flowing in it. Stepped leaders, 
dart leaders, and return strokes all propagate with different speeds 
because of the different physical processes occurring in them. Even 
within a given category of lightning phenomena the propagation speed 
can vary over a relatively wide range. This channel characteristic 
is important in that for propagation velocities much less than the 
speed of light, the aircraft response will begin long before the 
lightning current reaches the sensor point, because the electromagnetic 
fields from the current produce a polarization of the aircraft. 

3.1 Channel Modeling 

In order to calculate a realistic aircraft response, it 
is necessary to model the channel impedance and propagation velocity 
properly. The linear finite difference code used in this study 
achieves that by representing the lightning channel as a thin wire [4]. 
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In this approximation the channel is required to be much smaller than 
the cell size which for the FI 06 code is one-half meter. The thin 
wire formalism characterizes the channel impedance and propagation 
velocity through a per unit length inductance and capacitance. The 
following relationships hold: 



V = FL C l" 1 / 2 
v p LL jrt J 


(3.1) 


where 1 Q = characteristic impedance, (ft), 

Vp = propagation velocity (m/sec), 

= per unit length inductance (h/m), 

= per unit length capacitance (f/m). 

It should be noted at this point the per unit length inductance and 
capacitance in Equation (3.1) are not to be considered the real physical 
values for the lightning channel. The thin wire formalism is a mathe- 
matical construction which allows one to model a current path in the 
finite difference code which is smaller than the cell size. The per 
unit length inductance and capacitance used in the formalism are those 
appropriate for a coaxial cable with inner conductor radius a, and 
outer shield radius , where AS is the spatial grid dimension of 
the code. Thus C and L may be written, 

L = ^ 1t£ / s - n (§f) 


L 



„ / AS, 

tn ( H ) 


(3.2) 


In addition to specifying Z Q and V p , the thin wire formalism allows 
the inclusion of a resistance per unit length along the channel. In 
principle all of these quantities can be time varying, and some cases 
were run in which the resistance was allowed to vary. However, the 
majority of the cases reported in this chapter did not allow any 
time variation. By requiring time invariance of the channel, the 
problem of the aircraft response to a given excitation became amenable 
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to Fourier analysis. This greatly simplified the methodology used in 
finding a source to produce a given aircraft response. 

In addition to the thin wire formalism, two other methods 
were used to model the lightning channel. The first is simply a variation 
of the thin wire method, in which the exit channel is changed from a 
thin wire to a line of zeroed electric fields, corresponding to a 
perfectly conducting exit channel of dimensions one-half meter by 
one-half meter. This was done to see if the characteristics of the 
exit channel were important in the aircraft response. Slight differences 
were seen between the perfectly conducting and thin wire channels as 
can be seen in the data to be presented later. 

The second alternate method of modeling the channel was to 
simply force the injected current to be a particular value regardless 
of the channel characteristics. This method was used in the early 
part of the channel studies and corresponds to a very high impedance 
channel for which the current is unaffected by the presence of the 
aircraft. This is unlikely to be true physically and is unable to 
predict the correct relationship between time correlated response 
measurements on the aircraft. Results for this method will therefore 
be omitted from this report. 

3.2 Data Interpretation Based on a Single Measurement 

To illustrate how the channel study was done, one set of 
parameters will be analyzed in detail, with each step in the analysis 
explained. The entry channel for the study always ran from the outer 
boundary of the finite difference code to the nose of the FI 06 in a 
straight line. The exit channel ran directly to the boundary from the 
exit point. The source used to drive the code was a transparent current 
source located at the boundary of the code at the end of the entry 
channel. By "transparent" source is meant one that allows reflections 
coming back along the lightning channel from the F106 to pass through 
without reflection. 

The case to be discussed had inductance per unit length of 
3.33 x 10 7 h/m, capacitance per unit length of 3.33 x 10 -11 f/m, and 
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resistance per unit length of 50ft/m. The exit and entry channels both 
used the same parameters. The first step in analyzing the channel was 
to use a standard sine-squared pulse shape (shown in Figure 3.1) as a 
current source. This source was chosen so as to have significant 
frequency content at all relevant frequencies. The finite difference 
code was run with this source to determine the response of the FI 06. 

The forward D response is shown in Figure 3.2. Then the waveforms 
of Figures 3.1 and 3.2 were Fourier transformed to obtain the frequency 
content. The results are shown in Figures 3.3 and 3.4, respectively. 

Next the frequency transform of the FI 06 response was divided by the 
frequency transform of the source current to obtain the transfer 
function between source and response. The result is shown in Figure 3.5. 
This transfer function is a source independent quantity for the linear 
problem here considered. That is, regardless of the source current 
waveform, the Fourier transform of the response divided by the Fourier 
transform of the source results in the same transfer function. This 
knowledge was then used in conjunction with the measured data. The 
measured D response of flight 80-018 (Figure 3.6) was digitized and 
Fourier transformed, the result of which is shown in Figure 3.7. Then 
this transform was divided by the transfer function to give the Fourier 
transform of the current source needed to produce the response of 
flight 80-018. The transform of that source is shown in Figure 3.8. 

The frequency representation of that source current was then transformed 
back into the time domain, the results of which are shown in Figures 
3.9 and 3.10 on two different time scales. The derived time domain 
source current was then used to drive the finite difference code. This 
was done for two reasons. The first was to check that the source actually 
reproduced the measured waveform; that is, that there were no major 
errors in deriving the source. The second reason was to determine the 
actual injected current, which in general can be quite different from 
the source current. The forward D response of the aircraft is shown in 
Figure 3.11 and the injected current in Figure 3.12. Note that Figure 
3.11 does match the measured D quite closely, as it must if there were no 
errors in the analysis. The last step in the analysis is to use the 
source current in the finite difference code once again, but this time 
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Figure 3.2 Forward D-dot Response of F106 to Standard 
Current Pulse of Figure 3.1 

(L = 3.33 X 10' 7 h/m 

= 3.33 X 10" 11 f/m 

R fc = 50 ft/m) 
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Figure 3.4 Fourier Transform of D-dot Response of Figure 3. 
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Figure 3.5 Source Independent Transfer Function Derived as 
Explained in the Text 
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Figure 3.6 D-dot Response of Flight 80-018 
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Figure 3.7 Fourier Transform of D-dot Response of Flight 80-018 
(Flatched Area Indicates Closely SDaced Resonances) 
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Figure 3.8 Fourier Transform of Current Source Needed to Give 
Forward D-dot Response of Flight 80-018 
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Figure 3. 


9 Time Domain Representation of Current Source Needed 
to Give Forward D-dot Response of Flight 80-018 
(to 500 ns) 
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Figure 3.10 Time Domain Representation of Current Source Needed 
to Give Forward D-dot Response of Flight 80-018 
(to 25 usee) 
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Figure 3. 


Time (ns) 


1 Forward D-dot Response of F106 to Current Source oi 
Figure 3.9 






without the F106 present. That is, the lightning channel ran completely 
across the finite difference problem space. The current at the former posi 
tion of the nose of the FI 06 was monitored, and it is shown in Figure 3.13. 
In principle, this current should be an indication of the lightning environ 
ment. It is the current flowing in the 1 ightning channel in the absence of 
the aircraft. Figures 3.14 and 3.15 are the same as Figures 3.12 and 3.13, 
except on a larger time scale to illustrate peak values. 

Figures 3.16 - 3.27 show the injected current at the nose of 
the FI 06 and the current at that position without the F106 present for 
a variety of channel parameters as indicated on the figures. The currents 
are all appropriate to produce the measured D-dot response of Figure 
3.6. In all cases the exit channel ran from the base of the tail to the 
boundary. There are several cases in which the exit channel from the 
aircraft is listed as being a perfectly conducting wire. This does not 
affect the run that was made with the F 1 06 removed. The channel for that 
run had no perfectly conducting section and was a uniform thin wire 
across the entire problem space. 

Some things should be noted about Figures 3.14 - 3.27. First, 
all of the currents injected at the nose of the F106 are very similar, 
both in character and amplitude. This is a reflection of the fact that 
the transfer function between the current injected at the nose and 
the D-dot response is dependent almost entirely on the geometry of the 
aircraft, and only incidentally on the channel parameters. However, the 
current flowing at the same positions without the FI 06 present is 
strongly dependent on channel parameters, as expected. The amplitude of 
the injected current which produces the measured D-dot should also be 
noticed. This amplitude for a given D-dot will in general be a function 
of two variables, the peak amplitude of the D-dot record and the length 
of the record in time. That is, large injected currents must result in 
either large D-dot peaks or long time records. The record of flight 
80-018 is both largest in amplitude and longest in time of the measured 
D-dot responses. Hence it is reasonable to infer that the injected 
current calculated for that record is the largest seen so far in the 
measured D-dot's, so all lightning strikes which produced D-dot records 
must have had peak currents of less than 10,000 amperes. In fact, since 
most of the D-dot records are significantly smaller both in peak amplitude 


40 



Current (A) 


ORiGj-’fpv; - tj 

OF POOR QUALITY 



Figure 3.13 Current at the Position of the Injection Point for 
the Case in Which the FI 06 Was Not Present 
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Figure 3.14 Actual Injected Current at Nose of F106 to Obtain 
D-dot Response of Figure 3.12 (to 6 ys) 
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Figure 3.16 
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Figure 3.18 Current Injected Into F106 (L^ = 3.33 X 10”^ h/m 

C = 3.33 X 10' 11 f/m 

l 

= 0 n/m 

Exit Channel a Perfectly 
Conducting Wire) 
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Figure 3.19 Current at Injection Point Without F106 

(L = 3.33 X 10' 7 h/m 
& _n 

C = 3.33 X 10 f/m 
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= 0 a/m 

Exit Channel a Perfectly 
Conducting Wire) 
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Figure 3.20 Current Injected Into FI 06 (L^ = 1.67 X 10" 7 h/m 

C^ = 6.67 X 10' 11 f/m 
= 50 n/m 

Exit Channel a Perfectly 
Conducting Wire) 
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Figure 3.21 Current at Injection Point Without F106 

(L = 1.67 X 10“ 7 h/m 

C = 6.67 X 10' 11 f/m 

i 

= 50 s/m 

Exit Channel a Perfectly 
Conducting Wire) 




Current (A) 



Figure 3.22 Current Injected Into FI 06 
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Figure 3.23 Current at Injection Point Without F106 

(L = 5 X 10' 7 h/m 

C = 2 X 10~ 10 f/m 

R = 50 a/m 

Exit Channel a Perfectly 
Conducting Wire) 
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Figure 3.24 
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Figure 3.25 Current at Injection Point Without F106 

(L = 1.67 X 10‘ 7 h/m 
C = 6.67 X 10 11 f/m 
= 0 fi/m 

Exit Channel a Perfectly 
Conducting Wire) 
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Figure 3.27 Current at Injection Point Without F106 

(L = 5 X 10‘ 7 h/m 
* _i n 

C = 2 X 10 lu f/m 

= 0 O/m 

Exit Channel a Perfectly 
Conducting Wire) 




and duration, it is probable that a majority of the recorded strikes 
so far have had current amplitudes of less than a thousand amperes. 

3.3 Time Correlated Data Analysis 

Near the end of the contract period work began on the time 
correlated D-dot and B-dot measurements taken on the F106 aircraft 
during the summer of 1982. These records have similar features which 
allow one to recognize that both were produced by a single lightning 
event. There is no information, however, on which record led the 
other in time; that is, whether the B-dot or D-dot event occurred 
first. The methodology for analyzing these records is basically the 
same as for single records. A channel is chosen and a source deter- 
mined as previously described for a single record, either B-dot or 

D-dot. Then the analysis is done for the second of the time correlated 
measurements and a second source determined. In general the two 
sources will be different. This indicates that something about the 
chosen lightning channel is incorrect. By varying the channel para- 
meters and orientation it should be possible to bring the two calculated 
sources into agreement, so that the one remaining source will produce 
both of the measured records. At this point, though the calculated 
source is not mathematically unique, it is hoped that it can be re- 
garded as physically unique, in the sense that if a third response on 
the aircraft were measured, it would automatically be predicted by the 
calculated source. If this were not the case, the process of varying 
the channel characteristics would need to be redone until one source 
produced all three responses. 

The analysis of the time correlated measurements has 
concentrated on the records from flight 82-039 (Figures 3.28 and 3.29). 
because of their relatively large dynamic range and short record 
length. These factors contribute to accuracy and shorter computer runs 
respectively. Figures 3.30 - 3.37 give an example of the methodology 
used in dealing with the time correlated data. A channel was chosen 
and a source derived from the measured D-dot record. This was then 
used to drive the finite difference code. The resulting injected 
current is shown in Figure 3.30 and its time derivative in Figure 3.31. 
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Figure 3.28 D-dot Response of Flight 82-039 
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Figure 3.29 B-dot Response of Flight 82-039 
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Figure 3.31 Time Derivative of Injected Current for D-dot Source 
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Figure 3.32 Calculated D-dot Response for D-dot Source 
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Figure 3.36 Calculated B-dot Response for B-dot Source 
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Figure 3.37 Calculated D-dot Response for B-dot Source 
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The calculated D-dot response is given in Figure 3.32. Note that it 
corresponds quite well to the measured D-dot response. In Figure 3.33 
is shown the calculated B-dot response, which is clearly quite different 
from the measured response. Hence the chosen channel must be in error 
and needs adjusting. The reverse process was also done, in which a 
source was calculated for the B-dot record and used to drive the finite 
difference code. The injected current and its time derivative for this 
case are shown in Figures 3.34 and 3.35 respectively. Figure 3.36 gives 
the calculated B-dot response which again closely matches the measured 
response. In Figure 3.37 is shown the predicted D-dot response. The 
obvious difference between measured and calculated response verifies 
that the channel parameters chosen were in error. 

It should be pointed out that another possibility for error 
exists in addition to the possible error in the channel. That possibility 
is that the correlated responses were produced by a nonlinear event, 
in which case the model itself is in error and the analysis will be 
unable to derive a correct source. For this possibility a more sophis- 
ticated and yet undeveloped methodology will be needed. 

3.4 Time Varying Channels 

As mentioned earlier, some rudimentary work has been done for 
the case in which the channel resistance varied in time. An example of 
the results is shown in Figure 3.38 which is the calculated D-dot response 
for a step function current of amplitude 1000 amps and rise time 44 nano- 
seconds. The rear channel resistance per unit length varied as, 

r £ = — U -~- fl/m, I > 447 A 

= 500 n/m , I <_ 447 A. 

In addition the resistance was forced to be a monotonical ly decreasing 
function. That is, if the current I decreased, stayed at its minimum 
value. This models in a crude way the breakdown of the air as a lightning 
current flows through it. Since the model is not able to be analyzed by 
Fourier methods, little effort has been expended on it. More interesting 
results would be obtained by implementing real time varying channel 
models from the literature into the finite difference code. 
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CHAPTER 4 


NONLINEAR AIR BREAKDOWN MODELING AND RESULTS 
4.0 Introduction 

The interaction of lightning with an aircraft in flight is 
a nonlinear event. The initial attachment to the aircraft by the 
stepped leader is of necessity nonlinear, as is the development of 
an exit channel after the aircraft has charged to a sufficiently 
high level. Other aspects of the interaction, possible return 
strokes, K-changes, or dart leaders, also involve nonlinear processes, 
although it may be that the nonlinear effects are less important 
for these, because a conducting channel has been established previously. 
To accurately model these events, especially the initial leader, it 
is necessary to include electrical corona formation and air breakdown. 

It may also be necessary to take into account the effects of streamers. 

It is not difficult to understand in an elementary fashion 
what occurs as a lightning leader approaches an aircraft. The leader 
tip is a region of high charge density with correspondingly high 
electric field levels near it. As the leader approaches the aircraft 
the field intensity is largest in the direction of the aircraft. 

Free electrons in the air are accelerated in this high field until 
they collide with a neutral atom or molecule. If the electron's kinetic 
energy is large enough at the time of the collision, the neutral 
particle can have an electron separated from it, producing a second 
free electron and a positive ion. The free electrons are then again 
accelerated by the field, possibly suffering more collisions and 
producing more free electrons and ions. If the rate of production 
of free electrons is larger than the rate of loss (by recombination, 
and attachment to form negative ions), an electron "avalanche" 
occurs, in which sufficient numbers of electrons and ions are 
produced to substantially alter the electrical conductivity of 
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the air. The electrons and positive ions move in opposite directions 
under the influence of the electric field, constituting a current in 
the direction of the field. This has the net effect of neutralizing 
the electric field at the position of the leader tip, but enhancing 
the field at a point nearer the aircraft. That point now has a 
large charge density and is the new tip of the leader. 

The preceding description contains only what are commonly 
known as "primary" effects. There also exist "secondary" effects 
which may be important in lightning. 

The most significant of these is the photoionization process 
which allows for streamer propagation. A streamer is a luminous pulse 
of ionization which extends from the high electric field region into 
regions of lower field. The tip of the streamer is believed to consist 
of positively charged ions with a high enough density that breakdown 
or near-breakdown electric fields exist in some region around the tip. 
Also, some of these ions (which were created through collisions 
between electrons and neutral gas molecules) are likely to be in 
excited states which can decay and emit a photon. The energetics of 
the situation are such that this photon can ionize a gas molecule in 
advance of the streamer tip. The electron produced in the ionization 
can then avalanche in the field around the streamer tip, and the 
avalanche electrons drift back toward the tip leaving behind the 
positive ions formed in the avalanche. These positive ions become 
the new tip of the streamer. In this way a streamer can advance from 
regions of high field to regions where the field may be very low. 

4.1 Basics of the Nonlinear Model 

The model developed here to account for electrical corona 
is a variation of that commonly in use in nuclear electromagnetic 
pulse (NEMP) work for the calculation of air conductivity. The NEMP 
model solves for the air conductivity by calculating the densities 
of positive ions, negative ions, and electrons as a function of space 
and time through the use of detailed balancing. Physical processes 
included are electron avalanching, electron attachment to neutral 
molecules to form negative ions, electron-positive ion recombination. 
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and negative-positive ion recombination. 

The air conductivity for this model is a nonlinear function 
of the total electric field. It is given by 

a = q(n g y e + ( n _ + n +) • 

Here q is the charge on, one electron 1.6x10 19 coulombs, 

n is the number density of secondary electrons [m 3 ], 
e 

n and n + are the number densities of negative and 
positive ions [m~ 3 ], and 

u e and p.j are the electron and ion mobilities in 
[m 2 /(volt sec)]. 

The electron and ion densities are computed from the ambient 
ionization rate Q(t) : 

3n 

+ [3n + + a fi - G] n g = Q(t) 

3n 

if + n - = a e n e 


Here a is the electron attachment rate (sec 1 ), 
e 

G is the electron avalanche rate (sec *), 

0 is the electron-ion recombination coefficient (m 3 *sec 1 ) , and 
6 is the negative-positive ion recombination coefficient 
(m 3 -sec -1 ) . 

The coefficients are defined in Table 4.1. 
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TABLE 4.1 Air Chemistry Coefficent Formulas [5] 
Calculation of Erel : 


Erel = f- /(1+2.457P 0,834 ) for < 0. 07853(1+2. 457P U ’ 834 ) 


0.834, 


Erel = — -1 . 195P 0 ' 834 
p r, 


for — > 3.015+1 . 195P 0 * 834 
p r 


— +( — 
p r \ 


6884P 

2 


0.834' 


0.6884P 

2 


0.834 


for al 1 other — 
p r 


Where P is the percent water vapor and p^ is relative air 
density. Note: E is in esu, where Eesu = Emks/3xl0 4 . 

Calculation of Electron Attachment Rate a : 

e 

a e * (a 3 (l+0.344P)+a 2 ) 


a 2 = 1.22xl0 8 p r e' 21 ’ 15/Erel 

as = P r 2 (6.2xl0 7 +8.xl0 10 Erel 2 )/(l+10 3 Erel 2 (Erel (l+0.03Erel 2 ) ) ^ 3 ) 


Calculation of Avalanche Rate, G : 

G = 5.7xl0 8 p r YV(l+0.3Y 2 ’ 5 ); Y = 

Calculation of Electron-Ion Recombination Coefficient , and Ion-ion 
Neutralization Coefficient, : 

6 = 2x10 13 + p r 2.1xl0 -12 (m 3 /sec) 

B = 2x10” 1 3 + 2.8xl0” 12 (P) 1/3 (m 3 /sec) 

Calculation of Electron Mobility, 

100u 

“e ’ WTO '• R ■ 1 .55+210/(1+1 1 .8Erel+7.2Erel 3 ) - ^ s S ^ ter 

* (((16. 8+Erel)/(0. 63+26. 7£rel))°' 6 }/(3.xo r ) 

Calculation of Ion Mobility, u. : 

2.5X10” 1 * m/sec 
u i ” p p volts/m 
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Some of the weaknesses of the model are evident from the 
equations. The rate equations involving terms have no provision 
for motion of the charges, so electrons and ions must remain at 
the location where they were formed. This may be an adequate approxi- 
mation for very fast events such as those seen in NEMP work, but is 
not adequate for the relatively slow processes occurring in 1 ightning. 

Also the equations require local neutrality by specifying that the 
density of positive ions be equal to the sum of the densities of 
electrons and negative ions. This certainly is not true in the case 
of lightning. Finally, the physics of streamers is not included. 

In order to account for these difficulties, the equations 
have been modified in the following way. The rate equations involving 
terms and the local neutrality condition n + = n g + n_ have been 
replaced by a continuity equation for each species separately. These 
equations are shown below. 

an 

+ V • (n e v e ) + [Bn + + « e - G] n e = Q(t) 

an . , 

— rP- + v . (n_ v_) + 6n + n _ " a e n e (4.1) 

Sn + 

+ v • (n + v + ) + en e n + + 6n_n + = Q(t) + Gn g 

Note here that the source terms remain formally unchanged from the 
previous equations, and that the convective derivative has replaced 
the previous simple partial derivative, allowing for charge motion. 

The reason that the source terms are said to be formally unchanged is 
that the addition of streamer physics to the model changes the value 
of the avalanche rate G. Also a current of positive charge has been 
included in the convective derivative term to account for the effective 
motion of streamers. It should also be noted that the velocities 
v , v + , v in Equations (4.1) are calculated as the product of 
mobility and electric field, v s =+y $ E. The exception to this is the 
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case of v + when streamer propagation is included. 

4.2 Streamers in the Finite Difference Code 

The addition of photoionization and streamers to a finite 
difference code is a nontrivial matter, because the grid size is 
normally by necessity much larger than the radius of a streamer tip. 
Hence the code cannot sense the presence of streamers directly, because 
the effects of averaging over a spatial cell mask the streamer tip 
field. This makes it necessary to use special techniques and some 
assumptions in order to determine whether streamers exist in a given 
cell. The assumptions made are listed below. 

1) All net positive charge above a certain threshold is 
in the form of streamers. 

2) All streamer tips are in the form of spheres with an 
experimentally determined average radius and charge 
density. 

3) All streamers move with a uniform speed in the direction 
of the maximum electric field at the tip. 

None of these simplifying assumptions corresponds exactly to 
the physical situation, but they should be reasonable in an average 
sense. Since the finite difference code is dealing with quantities 
averaged over a spatial cell, it should be possible to get an accurate 
overall response from average streamer properties. 

The nonlinear finite difference code makes use of the above 
three assumptions in the following ways. Assumption (1) allows the code 
to determine which cells have streamers in them. These cells are then 
treated differently than cells without streamers. Assumption (2) allows 
the code to calculate the net positive charge density threshold and 
the number of streamer tips in a cell. That is, the experimentally 
determined radius and charge density of a streamer tip specifies the 
total positive charge in the tip. The net positive charge in a cell 
(volume of cell x(n + - n g - n_)) must be larger than this total positive 
charge in order for the cell to contain a streamer. If this is true 
then the total number of streamers in the cell is the total net positive 
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charge divided by the total charge in the standard streamer tip. 
Assumption (2) also allows the calculation of the electron avalanche 
occurring around the streamer tip. This is covered in more detail in 
the next section. Assumption (3) allows the code to model streamer 
motion by defining a current of positive charge in cells containing 
streamers. This current is proportional to the net positive charge 
density in the cell and also to the uniform speed chosen for the 
streamer. This uniform speed assumption is a weakness of the model, 
since streamers are known to move slowly as they form and then to 
speed up as they mature. 

4.3 Electron Avalanche Around a Streamer Tip 

In order to determine the effect a streamer has on the air 
conductivity it is necessary to calculate the magnitude of the electron 
avalanche which occurs around the streamer tip. To do this the streamer 
tip will be modeled as a uniform sphere of positive charge as shown 
in Figure 4.1. The charge density and radius will be assumed to be 
experimentally determined quantities. It will also be assumed that 
all avalanching occurs outside the streamer tip (i.e. r > r Q ). 



n i nr.21 -3 

n ps 1 x 1 0 m 

r Q « 2 x 10~ 5 m 


Uniform Positive Charge Number Density 

Radius of Streamer Tip 

Radius at Which Nominal Breakdown 
Electric Field Occurs 


Figure 4.1 Streamer Tip Geometry 
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The number of electrons being produced per unit volume, per 

unit time, around the streamer is G(r) n P (r), where G(r) is the avalanche 

P e 

rate, given in Table 4.1, and n g (r) is the number of photoelectrons 
per unit volume. The latter quantity is given [6] as 


n e P ( r ) = N s 




(4.2) 


where N $ = number of positive ions in the streamer tip, 
p = gas pressure, 
p = quenching pressure, and 

y _ number of photoelectrons generated/steradian/cm/Torr 
* number of ions formed in avalanching 


G(r) is given in Table 4.1 as a function of electric field. Because the 
field variation around the model streamer tip is known, this can easily 
be transformed to a function of radial distance. The result is; 

o (rv/r) 10 1 

G(r) = 5.7 x 10% —2 sec 1 (4.3) 

r l+.3(r b /r) 5 

The total number of electrons produced per unit time in the streamer 
avalanche is then just the integral of G(r) n e ^(r) over the volume 
surrounding the streamer tip. 


Hj streamer 



G(r) n g P(r) r^ dr. 


(4.4) 


The factor of one-half is included because the streamer tip is expected 
to be connected to a conducting channel, so high electric field and 
therefore the electron avalanche are confined to the half space in 
advance of the streamer. The number of electrons per unit volume per 
unit time produced in the finite difference cell by the single streamer 
is then, 


n 


e 


streamer 


streamer 

% 

'cel 1 


N s PP q t 
^cel 1 


1.14 x 



K/r) 


10 


r b 5 

l + .3(^) 


— dr, 
(4.5) 


where Ar is the approximate radius of the cell. The replacement 
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of the infinite upper bound in Equation (4.4) with Ar is a good 
approximation as long as Ar»r b> This is true in most cases. The 
integral in Equation (4.5) is easily transformed to give, 


streamer 


9 N s pp o t 
1.14 X 10* TT — - — — T ~~ - ~T 

V cell (p+p q ) 


J 


r b /r o 

r b /Ar 


u 


1 + .3 lT 


du. 


(4.6) 


The remaining integral can be performed numerically since r^, r , and 
Ar are known. 


Equation (4.6) gives the number of electrons per unit volume 
and time produced in a cell due to the presence of a single streamer. 
Assuming no interaction between streamers, the number of electrons 
per unit volume and time due to all streamers is. 


n all streamers = ^ n streamers 


(4.7) 


where N is the total number of streamers in a cell. To calculate N 

it is assumed that any net positive charge in a cell above a certain 

threshold is in the form of streamers. The threshold is the positive 

charge contained in a single streamer tip. So for the finite difference 

code to classify a cell as containing streamers, it must be the case 
C6 1 1 4 *5 

that ^net positive > 3" "o "o^’ follows directly from this that the 
total number of streamers in a given cell can be written as 


K - n e - n.) V „ 43 

- for („ + - „ e - nj V ce)) > j ,r n Q 


4/3 irr n 
o o 


for (n, - n_ - n ) V__ < i irr ^ n 

+ e cell 3 o o 


(4.8) 


^trp^mpr 

The quantity N n g then is the number of electrons per unit 
volume and time generated in a cell by all the streamer tips which 
are present. This term is then added to the right sides of the 


at 


and 


an, 

at 


equations in (4.1) as a source term. 


4.4 Streamer Motion 


The physical motion of a streamer tip is incorporated into the 
finite difference code as a current of positive charge. It couples into 
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the rate equations for the charge densities through the convective 
derivative term. The magnitude of this current is determined by the 
experimentally observed average streamer velocity, 2 x 10 m/s [7]. 

The current has the effect in the code of depleting a cell of any net 
positive charge, which stops the streamer avalanche in that cell. 

4.5 Two Dimensional Results 

The nonlinear corona model just described has been implemented 
in a two dimensional finite difference code which assumes cylindrical 
symmetry. This was done in order to model the rod-plane gap experiment 
of Collins and Meek [7]. Their experiment consisted of a positively 
charged rod and a negative plane as shown in Figure 4.2. The voltage 
applied between the rod and plane was such that electrical corona 
formed around the end of the rod. Time domain measurements were made 
of the electric field at the tip of the rod and on the plane directly 
below the rod. The results for a 64 kilovolt applied voltage are shown 
in Figure 4.3. The dotted line indicates the geometric field, which 
is the field that would be seen in the absence of corona. 


.6 m 



Figure 4.2 Rod-plane Gap Geometry of Collins & Meek [7 ]. 
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Figure 4.3 Measured Electric Fields for Experiment 
of Collins and Meek [7]. 

The calculated fields from the two dimensional nonlinear code 
are shown in Figure 4.4, overlaid on the Collins and Meek data. The 
major differences are in the final level of the field at the plane, and 
in the behavior of the field at the rod tip during corona formation. 

The excess of the measured field over the predicted field at the plane 
indicates that the experimental corona is larger spatially than that 
calculated by the computer code. It is likely that inaccuracies in the 
streamer formalism (e.g. average size, average velocity, etc.) included 
in the code are the reason for this. It should be nbted, however, that 
without the streamer model in the code, the calculated field at the 
plane rose just slightly ('v 5%) above the geometric field. 

The calculated and experimental fields at the rod tip differ 
in the amount of decrease seen during corona formation. The experi- 
mental field drops sharply to about .5 MV/m and then rises slightly 
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Figure 4.4 Comparison of Experimental Data of Collins 
and Meek [7] With the Prediction of the 
Nonlinear Computer Code 

before beginning a long slow fall. The calculated field also 
drops sharply but only to about 1 MV/m. From that point it continues 
to drop slowly, exhibiting no rise as is seen experimentally. It is 
felt that the cause of this difference is in the choice of propagation 
velocity for the streamers. It is known experimentally that this 
velocity is slower for streamers which are forming, and faster for 
mature streamers. Computationally, then, the use of an average velocity 
removes the streamers from the rod tip area too soon. Hence the 
electron avalanche around streamer tips occurs for a shorter time near 
the rod than it should, and the electric field stays artificially high. 
Presumably this defect in the model could be removed by a more appropriate 
variation of streamer velocity. 

Although Collins and Meek did not report results for a negatively 
charged rod, the prediction of the nonlinear code is shown in Figure 4.5. 
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Figure 4.5 Predicted Electric Fields for the 
Case of a Negatively Charged Rod 

Note that the field at the rod is much the same as in the positive rod 
case, but the field at the plane is essentially the geometric field. 

In this case, then, the spatial extent of the corona around the rod tip 
is much smaller than for the positively charged rod. This effect is 
confirmed by experimental observation. 

4.6 Application to the Three Dimensional Finite Difference Code 

Two applications were made of the nonlinear corona model to the 
three dimensional finite difference code. The first consisted of a 
perfectly conducting bar in free space and an electric charge near one 
end. This was done in order to study the effects of water vapor content 
of the air and relative air density on the attachment process. The 
second application was to study the attachment of an electric charge to 
the F106 aircraft. In this case the charge was placed near an extremity 
of the aircraft (i.e. wing tip, nose, tail tip) where lightning usually 
attaches. The objective of this application was to study the response of 
the F106 to nonlinear attachment. 
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It was found that for the cell sizes appropriate for these two 
problems { meter), the treatment of streamers is inaccurate. This 
occurs because of the way in which the streamer threshold (Equation 4.8) 
for a cell is calculated. In a cell with a small volume the difference 
n + -n -n_ must be large in order for Equation 4.8 to flag the cell as 
one having streamers. However, for the three dimensional codes of this 
chapter the necessary size of n + -n -n_ is of the same order as the 
roundoff error of the computer. Therefore, streamer cells cannot be 
accurately flagged. This problem could be resolved by using a smaller 
cell size or performing the computations in double precision, both of 
which require extraordinary computational resources. Hence, in the three 
dimensional results, the streamer formalism was not used, but the new 
terms involving charge motion were included. 

4.7 Parameter Study 

A parameter study was conducted with a small version of the 
nonlinear code in order to determine the effects of water vapor and air 
density on the attachment process. The size of the problem space was 
reduced to 8 m x 8 m x 29 m in order to decrease the running time of 
the code. A rectangular bar five meters long and one square meter in 
cross section was placed in the center of the space. A line current 
was then introduced into the space in such a way that an electrical 
charge appeared at a point two meters from the end of the bar. This 
charge was then allowed to grow in magnitude (as the time integral 
of the line current). Eventually the electric field between the charge 
and the bar reached air breakdown level, and a conducting channel 
appeared. The line current then could flow directly onto the bar. 

The parameters which were varied were the water vapor content 
of the air and the relative air density. The water vapor content is 
defined as the percentage of water vapor molecules by number in the air. 

The amount of water vapor in the air was allowed to be either 0 % or 6%, 
and the relative air density was either 1 or .5. Zero percent water vapor 
with relative air density of one corresponds to dry air at sea level, while 
six percent water vapor and relative air density of one-half corresponds 
to a thunderstorm at high altitude (^5.2 km). The line current used to 
form the charge was in all cases the same, and was chosen to be a step func- 
tion with sine-squared leading edge of amplitude 1000 amperes and rise 
time of 100 nsec. The outputs from the code were normal electric fields 
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and tangential magnetic fields on the bar, representative of charges and 
currents on its surface. 

The results of the parameter study are shown in Figures 4.6 - 

4.13 in groups of three. The (a) in each group is the normal electric 
field at the end of the bar facing the charge; (b) is the normal electric 
field at the end of the bar away from the charge; and (c) is the current 
at the end of the bar facing the charge. Note that this current includes 
a component due to polarization of the bar, which is also present in the 
linear case. 

There are a few things about Figures 4.6 - 4.13 which may 
appear puzzling at first glance. First of all, in some cases the fields 
seem to break down at levels which are lower than nominal values. 

This is entirely due to the "graininess' 1 of the finite difference 
code. Corona and air breakdown are, of course, most likely to occur 
around pointed objects such as the nose of an aircraft. These are 
regions of high electric field which causes the necessary electron 
avalanche. However, these regions in general also have the largest 
spatial gradients of the electric field. Therefore this field in a 
finite difference code can change markedly between two adjacent cells. 

One of the cells may have a field which is significantly above 

nominal air breakdown level and the other significantly below. The air 
conductivity in the nonlinear finite difference code is caluclated 
at spatial positions between electric field points. Hence for corona 
to form in the code it is necessary that the average of the fields 
in two adjacent cells reach breakdown level. So in the parameter study 
a field breaking down at a lower level than one would normally expect 
simply means that the field next to it was sufficiently high to make 
the average field between the two of breakdown intensity. 

Another characteristic of the parameter study which may be 
puzzling is the behavior of the field at the end of the bar facing the 
charge. In the case when the charge is positive the electric field 
first exhibits a slow rise followed by a sharp rise and then a steep 
fall when nonlinear levels are reached. The case of the negative 
charge is similar but shows no sharp rise. This is again partly due 
to graininess in the code, but another factor is also involved. Figure 

4.14 will help to clarify the situation. In (a) is shown the case of the 
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Figure 4.6c 
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Figure 4.8a 
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Figure 4.8b 
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Figure 4.8c 
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Figure 4.9a 
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Figure 4.10b 
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Figure 4.10c 
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Figure 4.11a 
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Figure 4.11c 
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Figure 4.12b 
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Figure 4.13a 
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Figure 4.13b 
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OP’f'yTvA L Ph r ' c iS 

OF POOR QUALITY 


(a) 


(b) 


Electron Flow 



Down First 


> Electron Flow 




Down First 


Fiq 4 14 Geometry of Breakdown Process in 
the Finite Difference Code 

positive charge and in (b) is the negative charge case. In both cases the air 

conductivity rises in the same place, between the charge and the end of the bar. 
The difference is in the direction of electron flow in the two cases. 

In (a) flow is toward the charge, which will neutralize the electric 
field near the charge but enhance the field at the end of the bar. The 
positive ions in the breakdown region will flow in the opposite direction. 

This will tend to neutralize the field at the end of the bar, but since 
the mobility of these ions is so much less than the electron mobility 
the effect will be very small at first. The enhancement of the field at 
the end of the bar will cause the breakdown region to extend all the 
way to the bar’s surface. This is what causes the sharp drop in the 
field after the enhancement is seen. 

In (b) electron flow from the breakdown region is in the 
opposite direction, causing the field at the end of the bar to be 
neutralized first. If one were to monitor the fields on the charge 
side of the breakdown region, an initial enhancement would be seen 
there, followed by a sharp drop. 

The question must now be asked as to whether the preceding is 
correct physical behavior or not, and whether or not the discreteness of 
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the finite differencing and the necessary spatial averaging have 
destroyed the physics of the attachment process. In the microscopic 
sense the answer is yes, in that the attachment does not take place 
exactly as one would expect. Physically, because of the cell size in the 
code, the bar has a relatively blunt end and also the maximum fields should 
be seen very near the charge. These fields should break down first pro- 
ducing an enhancement in the field at the end of the bar, regardless 
of the sign of the charge. But one must remember that the finite 
difference code cannot model in detail phenomena on a scale smaller 
than the cell size, and the initial air breakdown around the charge 
is surely smaller than that. The code can only model the breakdown 
in an average sense. That is, details such as the exact location of 
the initial electron avalanche are lost, but the overall behavior 
should be approximately correct. In particular, fields far from the 
attachment point should be relatively insensitive to the fine details 
of the breakdown and their behavior should be more characteristic 
of the average that the finite difference code calculates. 

One last characteristic of Figures 4.6 - 4.13 which may be 
confusing is the drop in the field at the end of the bar when attach- 
ment occurs. Notice that where a positive charge is present the 
field drops to a relatively low level (^2x10^ V/m) and stays there. 

In the case of a negative charge, however, the field drops to zero. 

This again is because of spatial averaging in the code and the 
direction of electron flow out of the breakdown region, as shown in 
Figure 4.14. The electrons are much more effective than the ions in 
neutralizing a field, so the fields toward which the electrons flow 
will drop to zero, and those on the other side will not. Again the 
average behavior of the fields is what is important to the attachment 
study. 

Some general comments can be made about the results of the 
parameter study. As expected, in all cases breakdown of the air around 
the charge is accompanied by a large pulse of current flowing onto 
the bar. This current rapidly drops after the initial breakdown to a 
much lower steady level. The magnitude of the current pulse 
is relatively insensitive to the water content of the air, but 
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is quite sensitive to the relative air density. This is because a low 
air density allows breakdown at a lower electric field intensity. The 
attachment, therefore, occurs at an earlier time when less charge has 
built up at the end of the bar. Hence a smaller current pulse is seen. 

The same type of phenomena is seen when comparing the current 
pulses for the two cases of positive and negative charge. The negative 
charge case shows an earlier and smaller current pulse for equal water 
content and relative air density. This occurs because the negative 
charge is composed of electrons which are available to start the electron 
avalanche earlier than in the positive charge case. So again, less 
charge has accumulated at the end of the bar when attachment occurs, 
and a smaller current pulse is seen. 

In none of the runs done for the parameter study was break- 
down observed at the back of the bar. Observation of the normal 
fields along the length of the bar at late time (^1 \isec) showed 
that the charge was being distributed fairly evenly, with less than 
a factor of two difference in charge density between the center of 
the bar and the back end. This implies that the formation of a con- 
ducting channel at the back end of the bar is likely to occur several 
microseconds after attachment. It is also likely to happen more 
slowly than the initial attachment, because there is no enhancement 
of the electric field as is seen between the charge and the bar in 
attachment. This is important to the F106 study in that it would 
appear that the initial attachment of a lightning channel to the air- 
craft and the initial exit will be quite widely separated in time, 
particularly for the small currents that are commonly seen. Of course, 
it should be kept in mind that the FI 06 aircraft has much more pointed 
structures, such as wing and tail tips, than the bar in the parameter 
study. These allow for much larger electric fields locally, and may 
cause exit channels to form significantly earlier than would be pre- 
dicted from a study of the behavior of the bar. 

4.8 Application of the Nonlinear Model to the F106B Aircraft 

To determine the response of the F106 aircraft to a lightning 
attachment, the nonlinear corona model was applied to the finite 
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difference model of the F106. The linear version of the F 1 06 model 
was described in detail in [1]. Addition of corona and air 
conductivity necessitated one change to the basic finite difference 
code. This was a reduction in the problem space size by three cells in 
each of the three coordinates, corresponding to a three meter reduction 
in the dimension along the length of the aircraft, and one and one- 
half meters in the other two dimensions. The problem space size for the 
nonlinear code was, therefore, 16 meters in the wing to wing dimensions, 
28 meters along the fuselage, and 9-1/2 meters vertically. The re- 
duction was necessary in order to bring the computer run time of the 
code down to an acceptable level, because air conductivity calculations 
are costly in both storage and time. Even with the reduction the cal- 
culations took approximately 30 CPU seconds per time step on a Data 
General MV8000 with 1 megabyte of central memory. 

Attachments were done at four points on the aircraft for 
both positive and negative charges. These points were the right side 
of the fuselage near the nose, the tip of the tail, and the leading 
edges of the right and left wings near the tip (Figure 4.15). The 
attachment was forced to occur by introducing a line current into the 
problem space and terminating it at a point one meter from the expected 
attachment point. This resulted in a charge building up at the end 
of the line current until the electric field between the charge and 
the aircraft reached breakdown level. Then a conducting channel formed 
between the charge and the aircraft allowing the charge and the line 
current to flow directly onto the aircraft. The waveform for the line 
current used to form the charge was a step function of amplitude 
100 amperes and sine-squared leading edge with rise time 100 nano- 
seconds. It was found that this current led to attachment at approxi- 
mately 500 nanoseconds and that the attachment was slow enough that 
the one nanosecond time step of the code could resolve it. Larger 
currents resulted in earlier attachment with very fast field changes 
which the code was unable to track with a one nanosecond time step. 
Smaller currents took an unacceptably long time to produce break- 
down field levels. In addition, 100 amperes is thought to be the 
average value of current in a stepped leader [8,9], which the nonlinear 
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Figure 4.15 FI 06 Model Showing Nonlinear Attachment 
Points and Sensor Points 
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code is attempting to model. The air parameters used in all of the 
attachments were a water vapor content of 6 percent, corresponding to 
thunderstorm and rainy conditions, and a relative air density of .5, 
corresponding roughly to the altitude (5.2 km) at which the F106 has 
been struck by lightning. Monitor points on the model aircraft were E and 
D-dot underneath the fuselage near the nose, H and B-dot on the 
fuselage above the right wing, and the current and its time derivative 
at the expected attachment point. The magnetic field monitor point 
was positioned so as to sense currents flowing longitudinally along 
the fuselage. The D-dot and B-dot monitor points correspond to those 
which are measured experimentally when the F106 is struck by lightning, 
so these responses can be compared directly with the measured data. 

The reader may be curious as to why one of the attachments 
was done to the side of the nose rather than directly onto the nose 
of the aircraft. The explanation has to do with the noncubical cells 
used in the nonlinear code. Figure 4.16 shows the geometric situation 
for a charge directly in front of the nose of the aircraft. The field 
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Figure 4.16 Illustration of Finite Difference Cell Containing 
Charge Located Directly in Front of Aircraft Nose 
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labelled E y should be the largest field because of the presence of the 
nose of the aircraft. However, because the unit cell is twice as long 
in that direction as in the two perpendicular directions, the field 
labelled E z is effectively closer to the charge than is E y . Therefore 
the nonlinear code finds that E z is larger than E y and when breakdown 
occurs, the charge moves predominantly in the directions perpendicular 
to E y . Physically the charge is attempting to achieve a lower energy 
state by becoming spherical, or at least as spherical as possible in 
a Cartesian coordinate space. But this expansion perpendicular to E y 
delays attachment and eventually causes the attachment to be a 
diffuse one over all points on the nose of the aircraft. Of course, 
this is not seen physically, and so must be avoided. By forcing the 
attachment to occur on the side of the nose, the expansion does not 
take place, and a much more localized area is involved. The response 
of the aircraft should not be affected drastically by this shift in 
attachment point. It is also observed that attachment from the side 
is a physically realistic situation. 

The aircraft responses for the attachment calculations are 
shown in Figures 4.17 - 4.24. Only a window around the time of attach- 
ment is shown in each case to emphasize the fast field and derivative 
changes that are seen. In all cases the attachment of a negative charge 
is seen to be somewhat gentler than that of a positive charge, as 
evidenced by smaller peak fields and time derivatives. This is consistent 
with the results in the study of the bar reported in the previous 
section. There is little difference in the current injected onto the 
aircraft at the four attachment points. Peak currents range from 450 
to 540 amperes with rise times of the order of 30 nanoseconds. Peak 
I-dot ranges from 3 x 10 10 to 8 x 10 10 amperes/second. 

As expected, there is greater variation in fields recorded at 
the other monitor points. Peak values of D-dot are largest for nose 
attachment, and peak values of B-dot are largest for right wing attach- 
ment, although the corresponding H field for that attachment is not 
largest. 

An explanation is necessary as to the difference in responses 
for left and right wing attachment, which should be symmetric. 
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Time (ns) Time (ns) 

Charge - Positive Charge - Positive 

Attachment Point - Nose Attachment Point - Nose 

Figure 4.17e. Injected Current at Attachment Point Figure 4.17f. I-dot at Attachment Point 
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Attachment Point - Nose Attachment Point - Nose 

Figure 4.18a. Electric Field at Forward Sensor Point Figure 4.18b. D-dot at Forward Sensor Point 
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Figure 4.20a. Electric Field at Forward Sensor Point Figure 4.20b. D-dot at Forward Sensor Point 
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Figure 4.20e. Injected Current at Attachment Point 
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.21a. Electric Field at Forward Sensor Point 
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Charge - Negative Charge - Negative 

Attachment Point - Left Wing Attachment Point - Left Wing 

Figure 4.22e. Injected Current at Attachment Point Figure 4.22f. I-dot at Attachment Point 
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Time (ns) T ™ e ( ns ) 

Charge - Negative Charge - Negative 

Attachment Point - Right Wing Attachment Point - Right Wing 

Figure 4.24a. Electric Field at Forward Sensor Point Figure 4.24b. D-dot at Forward Sensor Point 
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Figure 4.24e. Injected Current at Attachment Point 





Since the H and B-dot monitor locations are on the right side of the 
fuselage, these responses were expected to differ. The E and D-dot 
monitor locations, however, are also not exactly on the line of symmetry 
of the aircraft, because of finite differencing requirements. The actual 
location of these points is a quarter meter to the left of this line, 
and therefore approximately a half meter closer to the left wing 
attachment point. If these points were on the line of symmetry, no 
difference in E and D-dot response would be seen between left and right 
wing attachment. 

The calculated nonlinear aircraft responses were also compared 
with the measured responses from 1982. The 1982 responses were used 
because of their time-coordinated character, which permits one to be 
sure that they were caused by the same lightning event. In par- 
ticular, consider the measured D-dot and B-dot responses from Flight 
82-039, shown in Figures 4.25 and 4.26, and the negative charge nose 
attachment of Figures 4.18b and 4.1 8d . Although the actual wave- 
forms are different, the amplitudes and general character are 

very similar. By general character is meant lifetime of the response 
and also the time from peak to peak. In addition the calculated in- 
jected current of about 450 amperes is consistent within a factor 
of two of the linear current needed to produce the measured response, 
as found in Chapter 3. None of the other attachment points is able 
to reproduce both the measured D-dot and B-dot amplitudes, so it may 
be reasonably concluded that the measured records of Figures 4.25 and 
4.26 came from a nose attachment. The observed differences in the 
detailed waveforms may be envisioned as resulting from a slightly 
different attachment point, lightning channel orientation, or line 
current waveform. 

At present comparison of the measured data to nonlinear 
calculations is in a preliminary state. There is clearly no 
obvious technique for determining a current source such as that pre- 
sented in Chapter 3 for the linear interaction. It may be necessary 
to begin by developing a data base of calculated nonlinear responses 
in order to recognize trends and tendencies in the measured data. 
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Figure 4.26 Measured B-dot of Flight 82-039 
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One encouraging note, however, is that on the basis of the small number 
of calculations that have been done so far it appears that linear and 
nonlinear analysis give approximately the same answer as far as the 
amplitude of the lightning current is concerned. If this could be estab- 
lished as true in general, analysis of the measured data would be 
greatly simplified. 
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CHAPTER 5 


DISCUSSION AND RECOMMENDATIONS 
5.0 Data Interpretation State of the Art 

The research reported in the first four chapters is the 
continuation of initial efforts which were reported in [1], The principal 
objective has been to develop a methodology whereby the lightning 
environment can be obtained from the measured data. A linear and time 
invariant approach based on a combination of Fourier transform and 
three dimensional finite difference techniques has been developed and 
demonstrated. This approach can obtain the lightning channel current 
in the absence of the aircraft for various channel characteristic 
impedances and resistive loading. When applied to a single measurement, 
however, the environment so obtained is not unique. 

An interesting finding from this approach concerns the 
effect of the channel impedance on the response. From an environment 
interpretation point of view, it was found that the current waveform 
injected onto the aircraft is not greatly affected by the channel 
parameters, but that the channel current in the absence of the aircraft 
is. 

Toward the end of this effort, time correlated 1982 B and 
D measurements were made available, and the approach was applied to 
one such data set. The results indicate that the initial selection 
of channel parameters (attach point, characteristic impedance, resistive 
loading, etc) did not reproduce both measurements simultaneously. There- 
fore more work needs to be done to understand simultaneous data. 

In order to help with this problem the concept of response 
ratio was introduced. Ratios were computed for some measured and 
computed responses. Results are encouraging, although the concept 
needs to be more fully developed and applied. 

It should be emphasized, however, that at this point there is 
no guarantee that linear analysis is sufficient to understand the 
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measured data. With this in mind, a nonlinear model of air breakdown 
based on air chemistry, photo-ionization, charge motion, and streamer 
theory has been developed and successfully compared with data published 
in the open literature. This model has been incorporated into the three- 
dimensional finite difference electromagnetic interaction model of 
the F106B aircraft. Attachment studies have been done for various 
parameters, and the results yield waveforms quite similar to, but not 
exactly the same as, those measured in flight. Indeed at the present 
time, it seems likely that the nonlinear attachment process could account 
for most of the in flight data obtained so far. More research and more 
simultaneous data is required to resolve the issue of whether or not 
the interaction can be successfully interpreted by linear analysis alone, 
or whether nonlinear modeling is required. 

It should be noted here that the emphasis thus far has been 
on methodology and tool development. This includes the nonlinear model, 
the linear method development, and the concept of the response ratio. 

A major effort in future research will be to apply the methods to the 
data, especially the 1982 time correlated measurements. 

5.1 Recommendations for Future Test Programs 

a. Increased Dynamic Range 

It has been shown in Chapter 2 how the digitization error 
limits the dynamic range of the data and thus its utility and accuracy 
in the data interpretation processes. Previous data has been obtained 
with a 6 bit recorder, and it is understood that an updated 8 bit 
version is planned. This will increase the available dynamic range 
by up to 12 dB, which should be a noticeable improvement. 

b. Several Time Correlated Measurements 

The 1982 data offers great hope in that B, D, and I records 
caused by the same lightning event were obtained. However, the precise 
timing of these records with respect to each other is not known. It is 
desirable that several truly simultaneous measurements be made at 
widely separated (in space) measurement points. The relative timing of 
the records should be accurate to at least 10 ns, so that it can be 
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determined which record began first. Knowledge such as this would be 
greatly helpful in identifying attach and exit points. 

c. Measurements of Fields 

The data obtained on the surface of the aircraft consists of 
D and B measurements. It is difficult to determine the initial condition 
for these variables, but this makes a difference in the interpretation 
process. It is useful to know, for example, if the decrease in electric 
field is from a high level towards zero, or if it is from zero to a 
negative value. Thus, knowledge of the field levels themselves would 
be of great help. 

d. Direct Strikes at Low Altitudes 

Most of the data obtained thus far has been for relatively 
small events at high altitudes. It is certainly desirable to obtain 
data at lower elevations so that return stroke events can be recorded. 


145 


REFERENCES 


1. Rudolph T.; and Perala, R.A., "Interpretation Methodology and 
Analysis of In-Flight Lightning Data," NASA CR-3590, October 1982. 

2. Pitts, Felix, L.; and Thomas, Mitchel E.: 1981 Direct Strike 
Lightning Data. NASA TM-83273, 1982. 

3. Turner, C.D.; and Trost, T.F.: "Laboratory Modeling and Analysis 
of Aircraft-Lightning Interaction," NASA CR-169455, August 1982. 

4. Merewether, D.E.; and Fisher, R.: "Finite Difference Solution 
of Maxwell's Equation for EMP Applications." Electro-Magnetic 
Applications, Inc. (EMA), P.O.Box 8482, Albuquerque, NM 87198, 
Defense Nuclear Agency, DNA 5301F, April 1980. 

5. Radasky, William, A.: "An Examination of the Adequacy of the 
Three Species Air Chemistry Treatment for the Prediction of Sur- 
face Burst EMP." Defense Nuclear Agency, DNA 3880T, Dec. 1975. 

6. Gal 1 imberti , I.: "The Mechanism of the Long Spark Formation," 
Journal de Physique, Collogue C7, supplement au No. 7, Tome 40, 
Juillet 1979. 

7. Collins, M.M.C.; and Meek, J.M. : "Measurement of Field Changes 
Preceding Impulse Breakdown of Rod-Plane Gaps," Proceedings of 
the Seventh International Conference on Phenomena in Ionized 
Bases, August 22-27, Gradevinska Kujiga Publishing House, 
Belgrade, 22-27 August 1965. 

8. Uman, M.A. ; and Krider, E.P.: "Lightning Environment Modeling," 
Volume I of "Atmospheric Electricity Hazards Analytical Model 
Development and Application." Electro Magnetic Applications 
EMA 081-R-21, and AFWAL-TR-81 -3084, August 1981. (Available from 
DTIC as AD A114 015.) 

9. Uman, M.A. ; and Krider, E.P.: "A Review of Natural Lightning: 
Experimental Data and Modeling," IEEE Transactions on Electro- 
magnetic Compatibility, Vol . EMC-24, No. 2, May 1982. 


146 



APPENDIX A 


ORIGINAL FACE ri 
OF POOR QUALITY 


COMPARISON OF EXPERIMENTAL ANO NUMERICAL 
RESULTS FOR THE INTERACTION OF A SCALE MODEL 
AIRCRAFT WITH A SIMULATED LIGHTNING CHANNEL* 

R. A. Perala and T. H. Rudolph 
Electro Magnetic Applications, Inc. 

P. 0. Box 26263 
Denver, Colorado 80226 

and 

T. F. Trost and C. D. Turner 
Department of Electrical Engineering 
Texas Tech University 
Lubbock, Texas 79409 


Summary 

Results are given which compare measured and 
computed responses for a scale model aircraft in a 
simulated attached lightning channel. A scale model 
of the NASA Langley F106 lightning research aircraft 
was suspended by a wire simulating the channel. A 
pulse was injected on the wire and subsequently inter- 
acted with the aircraft model. Sensors on the model 
recorded surface magnetic and electric field deriva- 
tives. Numerical simulations of this configuration 
were made with a three dimensional finite difference 
code for four entry/exit point configurations. Good 
agreement was obtained. 


Introduction 

One of the items of current interest in the 
lightning community is the knowledge of lightning's 
electromagnetic properties during a lightning/aircraft 
interaction event. This is of great importance because 
of the aircraft industry's move to new technologies 
such as digital fly by wire, advanced composite mate- 
rials, and extensive use of low level integrated cir- 
cuitry. Upset of digital circuits by lightning induced 
transients is therefore becoming a topic of increasing 
Interest. Because of this concern, the NASA Langley 
Research Center has been sponsoring a research program 
to investigate the electromagnetic characteristics 
of natural lightning at aircraft altitudes. The prin- 
cipal means of accomplishing this investigation has 
been their FI 06 thunderstorm research aircraft which 
has been Instrumented with electromagnetic sensors and 

recorders. * The aircraft is flown into thunderstorms 
with the intent of being struck by lightning. Data on 
the interaction of the aircraft with both attached 

and nearby strikes has been obtained. 3 ' 6 In order to 
properly interpret the data, that is, to identify the 
nature of the lightning that cause the aircraft re- 
sponse, two accompanying supportive parallel efforts 
are being conducted. The first approach is experimen- 
tal. A scale model of the F106 was constructed and 
suspended by wires which simulate the lightning chan- 
nel. A pulse was injected up the wire and aircraft 
surface magnetic and electric field derivative respon- 
ses were measured. The second approach involves numer- 
ical simulation of the scale model aircraft interaction 
by use of three dimensional finite difference solutions 
of Maxwell's equations. The results from the two ap- 
proaches are then compared in order to enhance the 
understanding of the interaction process. 


«»^ k ,^ 0rmed under NASA r,rant NAG 1-28 and contract 
NASI -1 6489 and subcontract 1-43U-2094 to Research 
Triangle Institute. 
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Experimental Approach 

An approximate scale model of the F-106B delta- 
wing aircraft has been constructed. 7 The fuselage con- 
sists of a 3 foot length of aluminum cylinder, 4 in- 
ches in diameter. Flat end caps are machined to fit in- 
to each end and are secured with screws. The wings and 
tail are cut from 1/16 inch brass sheet to scale with 
the aluminum cylinder. The overall scale of the model 
is 1/18.8 of the actual F-106B. 


The modeling apparatus is shown in Figure 1. 



Figure 1 Apparatus for Aircraft- 
Lightning Modeling 


A pulse of 0.75 ns duration is launched at the bottom 
of the lower wire. The magnetic (B) field measured six 
inches from the wire and 60 inches away from the 
model s nose is shown in Figure 2. To measure the re- 
sulting transient electromagnetic fields, time deri- 
vative sensors are used. Small B-dot (for longitudinal 
and transverse currents) and D-dot sensors have been 
placed on the model in the locations corresponding 
to their actual locations on the F-106B (on the fuse- 
lage over the starboard (axial current) and port wings 
(transverse current) and under the nose). The B-dot 
sensor is a loop made by bending 0.141 inch diameter 
semirigid coaxial cable into a 0.9 cm radius semi- 
circle and cutting a gap in the outer conductor. The 
D-dot sensor is a monopole using a 1.65 cm wire. 
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The pulse travels up the lower wire, over the 
model, and on up the outer conductor of the upper wire, 
which is actually 0.141 coax. The sensor outputs are 
carried on the inside of this cable. No measurable 
leakage to the inside of the cable has been detected. 
When the top of the apparatus is reached, a transition 
is made to 0.5 inch diameter cable to complete the run 
back to the sampling oscilloscope. There, the waveforms' 
are digitized and stored on magnetic disk. The equi- 
valent risetime of the sampling scope is 25 ps. Although 
Figure 1 shows the attachment points in the center of 
the end caps, results have been obtained for other con- 
figurations as discussed in the last section. 



Figure 2 Magnetic Field Measured at a 

Radial Distance of 6 Inches from 
the Wire and 60 Inches from the 
Model Airplane's Nose 

Numerical Approach 

The numerical approach involves the use of a 
three dimensional time domain finite difference solu- 

9-12 

tion of Maxwell's equationsin Cartesian coordinates. 

The relationship of the numerical model geometry to the 
actual shape is shown in Figure 3. 



The spatial resolution of the solution space is one 
meter in the long dimension of the aircraft and one- 
half meter in the other two dimensions. The temporal 
step size is one nanosecond. The actual size of the 
solution space was chosen so that there are at least 
eight cells between any point on the aircraft and the 
boundary. The boundary conditions used on the outer 
boundary are the radiation boundary conditions of 

Merewether. 12 Current was injected on the aircraft 
by specifying magnetic fields around the position of 
the input wire, and removed from the aircraft by zero- 
ing a line of electric fields from the plane to the 
boundary at the position of the exit wire. 

Figure 2 , from which the injected current was 
deduced, needs further explanation. The measurement was 
taken at a point five feet from the nose of the model, 
and six inches from the wire, so the B field in the 
plot is an indication of the current at that point in 
the wire,- which is not the same as the injected current. 
That this is true can be seen from the waveform of the 
B field, which shows a second peak approximately ten 
nanoseconds after the first peak. The ten nanosecond 
period corresponds to the travel time of the signal 
from the measuring point to the model and back, indi- 
cating that the second peak in Figure 2 is a reflection 
from the model travelling back along the wire. A re- 
flection occurs at the injection point of the model 
because of the mismatch of impedances there. The model 
presents a lower impedance to the current than the wire 
does, and hence, more current is injected onto the mo- 
del than is incident from the wire. This results in the 
reflected wave seen in the second peak of the plot. 

The actual injected current must then be determined 
from the sum of the two peaks, and not by the waveform 
of Figure 2 as it stands. For these studies the inject- 
ed current was deduced by shifting the maximum of the 
second pulse of Figure 2 to the time position of the 
first maximum and then summing point by point. 

A difficulty which arose in the analysis was in 
how to accurately determine the injected current from 
the magnetic field measurement. The most obvious way 
is to assume the simple expression B(t) = u 0 I (t )/2*r, 

where r is the perpendicular distance of the measure- 
ment point from the wire. But this is really a magneto- 
static assumption and at least requires that the pulse 
width be a much greater than the signal travel time 
from the wire to the measuring point. In the present 
case the pulse width is about .75 nanoseconds and the 
travel time is .5 nanoseconds, so the requirement is 
not satisfied. In order to solve the problem more accu- 
rately, an integral expression can be derived from the 

current in terms of the magnetic field. ^ This approach 
is complex, and it would have required a significant 
amount of time to unfold the true current from the 
measured field. It was determined that the cost in time 
was not justified and the decision was made to use the 
simple formula B(t) - u 0 I(t)/2Trr, even though it may 

not be as accurate as one would like. Further measure- 
ments of the magnetic field closer to the wire are 
planned, but are not yet available. 

Results and Conclusions 

Results were compared for the four aircraft 
entry - exit point configurations shown in Figure 4. 

A comparison of the amplitudes of the initial peak is 
given in Table 1. Overlays of measured and computed 
results are shown in Figures 5 and 6. 

The results show that in every case there is at 
least reasonable agreement between waveform shapes 
and in several cases the agreement is quite good. 
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Looking specifically at the measured waveforms 
of B l and D for the nose-tail and nose-wing cases, one 

can observe reflections of the incident pulse from the 
trailing edge of the wing and from the rear of the 
fuselage (and trailing edge of the tail). These occur 
about 2.5 and 3.5 ns, respectively, after the incident 
pulse in and about 5.0 and 5.2 ns after the pulse 

in D. As the rear wire is moved from the fuselage to 
the wing tip, the reflections from these two locations 
are altered, the reflection from the trailing edge of 
the wings changing from being larger to being smaller 
than the reflection from the fuselage. This is expect- 
ed because the wire carries away some of the current 
from the point to which it is attached, thus giving 
a smaller reflection. 



Nose-Tail 

Injection 


Nose-Wing Tail-Nose Wing-Wing 

Injection Injection Injection 


Figure 4 Aircraft Entry-Exit Point 
Conf i gurations 


The main sources of error in the measured fields 
include the uncertainties in amplitude and time-base 
calibration of the sampling oscilloscope, the uncertain- 
ties in the sensitivities of the D-dot and B-dot sensors 
and the bandwidth limitations of the sensors and co- 
axial cables. The sensitivities of the sensors have been 
measured in a biconical transmission line calibrator 
which was constructed for this purpose. Smaller errors 
are introduced by digitization (12 bit A/D converter), 

timing, and noise (S/N > 40 d8^). 

Errors in the numerical approach chiefly come 
from: 1) imperfect boundary conditions at the finite 
difference mesh outer boundary, 2) uncertainty in the 
exact knowledge of the input current waveform, 3) 
approximations associated with the ability to exactly 
model the aircraft shape with cell sizes 1m x .5mx ,5m. 
This also indicates some uncertainty in computing re- 
sponses at exactly the corresponding locations used 
for the measured data. 

The results are therefore felt to generally com- 
pare within experimental and numerical errors. The 
comparison therefore gives confidence that these 
methods can be applied to the F-106B data interpre- 
tation problem with the hope that valid conclusions 
about lightning interaction with aircraft can be made. 
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Table 1 Comparison of Amplitudes of Initial 
Peaks, and Ratio in dB 


ttst-nu 
feu lie jl 

560 393 9.3 

V 10 '* 1 ’ m 192 i.3 
us m j.; 

w | 

i(C/^.,) 7 2 5 , 1 7 

MID* 5 C/k 7 ) 2.2 2.S 1.5 


Olu Iuh * feu 

5*0 393 3.5 250 

227 ]« t,| 17t 

122 192 3.9 150 
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7.2 5.9 1.7 3.1 
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L -WQSt wmt-inw 

bn OL feu Km 
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77 9.1 1N0 101 ?.» 

165 .1 2*0 195 1.9 

1« 9.1 112 107 .1 

.97 11.9 2.5 7.0 1.9 

.1) 7.9 1.1 1.7 3.1 


The results for the calculated B^ and D waveforms 

show less change in the reflections between these two 
cases evidently indicating somewhat less sensitivity 
to the location of the rear wire. For B^, the reflect- 
ion from the rear of the fuselage is decreased a bit 
in the nose-tail case, and for D there is no change. 

The amplitude comparisons of Table 1 show an. 
average agreement of 3.6 dB. It is noted that the Bj 

and comparisons are usually the most divergent, 

which may indicate an unusual difficulty in the measure- 
ment or numerical technique at that location. The 
large differences of 11 dB are in two cases only. If 
these two 11 dB differences are not included, the aver- 
age value of the absolute difference between measured 
and predicted peak amplitudes is 2.9 dB. 

There are errors in both the measured and numeri- 
cal results, although it is generally difficult to 
assign exact values. The measurement error in amplitude 
is expected to be less than 2 dB, and one could easily 
account for similar values for the numerical technique. 
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APPENDIX B 


FOURIER TRANSFORM TECHNIQUE 

The convention adopted for the Fourier transform used in the 
analysis is the following: 


Time to frequency: 

F (ui) = 

/ 

— 00 

f (f) e‘ jujt, dt' 

(B-l) 

Frequency to time: 

f(t) - 

1 _ 

2 TT 

00 • 1 
f F(u' ) e Ja) t do,' 

(B-2) 


— oo 


For the time to frequency transform it is assumed that f(t) 

is nonzero for only a finite time, 0<t<x, and is represented in this 

interval as a sequence of uniformly spaced points; f Q , f -j , f^, ... , f^, 

where the uniform temporal spacing is At. It is also assumed that 

the sequence of points is chosen such that f = f N = 0. The function 

f ( t ) is constructed from the points f , f 1 , ..., f N , by connecting the 

points with straight lines. Therefore f(t) is a continuous function 

with discontinuous derivative at each of the points f , f , , .... f„. 

o 1 N 

Equation (B-l) can be integrated by parts to yield, 

F( W ) = - f 3 ( ^e- Ju)t ' |* +3I f] f(t') e- jwt ' dt\ ( B-3) 

But f ( x ) = f N = 0 and f(0) = f Q = 0, so, 

F <“> = -J ’l f ' (*■) •" M ' «'• (8-4) 

Breaking the interval 0<t<x into segments of length At, 

Hji f , (t .) ,-j»r dt , ( 

w n=0 t 


Now f (t ) is a constant for each interval and can be removed 
from the integral. 


where f 


n 


, N-l x n+l . . , 

FU) = -J- f / e' Jut 

o) n-U n 

f - f n 

= V] 

At 


dt' , 


(B-6) 
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The remaining integral can now be performed to give, 


, N-l -jwt 

F <“>- i A f n' <- fc > (e 


n+1 


-jcot 


-e 


n ) 


N-l 


•jut. 


•k n5o 

“ , N -’ -Jut 

F(u) - m (e‘ J -') n ? 0 - V e 
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Note that t n = nAt and e 
so that, defining x= e Jt0 , 


-juAt a cons tant at each frequency, 


F(w) = 


x-1 


N-l 


L \l " f ) X 11 . 

At n=0 n+1 


(f 


(B-8) 


The transform has thus been reduced to the evaluation of a 
polynomial at each frequency. Equation (B-8) is the form of the time to 
frequency transform which has been implemented for the analysis in this 
report. It should be noted that F(u) as represented in Equation (B-8) is 
exact and known at all frequencies. This has come about because of the 
definition of f(t) as a collection of straight line segments, which makes 
f(t) known at all times. In reality the value of f(t) is completely 
unknown except at the discrete points f Q , f -j , ..., f N - Hence, to avoid 
aliasing. Equation (B-8) should be trusted only in the range w <_ ir/At. 
This ensures that at least two of the discrete known points of f(t) will 
fall within the period 2 tt/oj. 

The frequency to time transform. Equation (B-2), is implemented 
in much the same fashion as the time to frequency transform. It is comp- 
licated by a pair of things. First, F(w) is in general nonzero as w -*• °°, 
so the integration in Equation (B-2) cannot in principle be truncated 
at some large However, in practice, the contribution of frequencies 
beyond wu to the integral is usually very small, so Equation (B-2) is 
approximated by, 

f(t) = / F(u)‘ ) ^ da) 1 f F(u) 1 ) du> . (B-9) 

-“ N 
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The second complication is that the discrete frequency points 
one is usually dealing with, F Q , F-j , F^, are not in general uniformly 

spaced. This means that the final form of the frequency to time transform 
will not be as simple as that of the time to frequency transform. 

To implement the frequency to time transform assume that the 
discrete points are known: F Q , F ] , F 2> F^, where F q = F(u> ), 

F 1 = F (“-|)> .... F^j = F(«o n ). Also choose large enough so that F^ % 0. 


The first step in the implementation of Equation (B-9) is to 
eliminate the integration over the negative frequencies. Note first that 
from Equation (B-l), F(-u) = F(u>), where the bar indicates complex conjugate. 


f(t) 
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k f -« e 

i7T 

-I 0 . , 


p <“’> eJW 1 p (“') e jVt 
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(B-10) 


But F (a)') e‘ Ja,,t + F(a.') e jtl) ' 1 = 2 Real {F( u ') e jVt }. 


Thereto re, 

, 1 "N 

f(t) = - Real { f Q F(w' ) e Jw t dui 1 }. ( 

From this point the development proceeds as in the time to 
frequency transform. 
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This is the form of the frequency to time transform that has 
been implemented for the analysis in this report. Note that Equation (B- 
is really an approximation because of the truncation of the integral 
in Equation (B-2) at 
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